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Abstract 

We discuss three related subjects well suited to graduate research. The first, Nonequilibrium 
molecular dynamics or "NEMD" , makes possible the simulation of atomistic systems driven by ex- 
ternal fields, subject to dynamic constraints, and thermostated so as to yield stationary nonequi- 
librium states. The second subject, Smooth Particle Applied Mechanics or "SPAM", provides 
a particle method, resembling molecular dynamics, but designed to solve continuum problems. 
The numerical work is simplified because the SPAM particles obey ordinary, rather than partial, 
differential equations. The interpolation method used with SPAM is a powerful interpretive tool 
converting point particle variables to twice-differentiable field variables. This interpolation method 
is vital to the study and understanding of the third research topic we discuss, strong Shockwaves in 
dense fluids. Such Shockwaves exhibit stationary far-from-equilibrium states obtained with purely 
reversible Hamiltonian mechanics. The SPAM interpolation method, applied to this molecular 
dynamics problem, clearly demonstrates both the tensor character of kinetic temperature and the 
time-delayed response of stress and heat flux to the strain rate and temperature gradients. The 
dynamic Lyapunov instability of the shockwave problem can be analyzed in a variety of ways, 
both with and without symmetry in time. These three subjects suggest many topics suitable for 
graduate research in nonlinear nonequilibrium problems. 

PACS numbers: 05.20.-y, 05.45.-a,05.70.Ln, 07.05.Tp, 44.10. +i 
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I. THERMODYNAMICS, STATISTICAL MECHANICS, AND NEMD 



A. Introduction and Goals 

Most interesting systems are nonequilibrium ones, with gradients in velocity, pressure, 
and temperature causing flows of mass, momentum, and energy. Systems with large gra- 
dients, so that nonlinear transport is involved, are the most challenging. The fundamental 
method for simulating such systems at the particle level is nonequilibrium molecular dynam- 
ics (NEMD). 1-3 Nonequilibrium molecular dynamics couples together Newtonian, Hamilto- 
nian, and Nose- Hoover mechanics with thermodynamics and continuum mechanics, with the 
help of Gibbs' statistical mechanics, and Maxwell and Boltzmann's kinetic theory. Impulsive 
hard-sphere collisions or continuous interactions can both be treated. 

NEMD necessarily includes microscopic representations of the macroscopic thermody- 
namic energy E, pressure and temperature tensors P and T, and heat-flux vector Q. 
The underlying microscopic-to-macroscopic connection is made by applying Boltzmann and 
Gibbs' statistical phase-space theories, generalized to include Green and Kubo's approach 
to the evaluation of transport coefficients, together with Nose's approach to introducing 
thermostats, ergostats, and barostats into particle motion equations. 

These temperature, energy, and pressure controls make it possible to simulate the be- 
havior of a wide variety of nonequilibrium flows with generalized mechanics. The nonequi- 
librium phase-space distributions which result are typically multifractal, as is illustrated 



here with a few examples taken from our website, [ http://williamhoover.info ]. These 
ideas are summarized in more detail in the books "Molecular Dynamics" , "Computational 
Statistical Mechanics" and "Time Reversibility, Computer Simulation, and Chaos". The 
one-particle "Galton Board" (with impulsive forces) and the "thermostated nonequilibrium 
oscillator" problem (with continuous forces) are simple enough for thorough phase-space 
analyses. Macroscopic problems, like the steady shockwave and Rayleigh-Benard flow, can 
be analyzed locally in phase space by computing local growth rates and nonlocal Lyapunov 
exponents. 

The main goal of all this computational work is "understanding" , developing simplifying 
pictures of manybody systems. The manybody systems themselves are primarily compu- 
tational entities, solutions of ordinary or partial differential equations for model systems. 
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Quantum mechanics and manybody forces are typically omitted, mostly for lack of com- 
pelling and realistic computer algorithms. There is an enduring gap between microscopic 
simulations and realworld engineering. The uncertainties in methods for predicting catas- 
trophic failures will continue to surprise us, no matter the complexity of the computer models 
we use to "understand" systems of interest. 

Number-dependence in atomistic simulations is typically small: 1/N for the thermody- 
namic properties of periodic iV-body systems, perhaps 1/yN or even 1/lniV in problems 
better treated with continuum mechanics. So far we have come to understand the equilibrium 
equation of state, the linear transport coefficients, the Lyapunov instability of manybody 
trajectories, and the irreversibility underlying the Second Law. Improved understanding of 
relatively-simple hydrodynamic flows, like the Rayleigh-Benard flow treated here, will fol- 
low from the special computational techniques developed to connect different length scales. 
Smooth Particle Applied Mechanics, "SPAM" has proved itself as not only a useful sim- 
ulation technique for continuum systems, but also as a powerful interpolation tool for all 
point-particle systems, as is illustrated here for the shockwave problem^-. 

B. Development of Molecular Dynamics at and Away from Equilibrium 

In the early days of expensive vacuum-tube computing the hardware and software were 
largely controlled by the Federal Government and located at the various weapons and energy 
laboratories at Argonne, Brookhaven, Livermore, Los Alamos and Oak Ridge. Fermi devel- 
oped molecular dynamics at the Los Alamos Laboratory in the summers of 1952-1953, discov- 
ering many of the interesting nonergodic recurrence features characterizing the low-energy 
behavior of one-dimensional anharmonic "Fermi-Pasta-Ulam chains" . The Los Alamos Re- 
port summarizing his work was prepared a few months after his death^— . At sufficiently 
low energies the anharmonic chains showed no tendency toward equilibration while (it was 
discovered much later that) at higher energies they did. Figure 1 shows time-averaged 
"mode energies" for a six-particle chain with two different initial conditions. In both cases 
the nearest-neighbor potential generates both linear and cubic forces: 

0(r) = (l/2)(r-l) 2 + (l/4)(r-l) 4 . 
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FIG. 1: The time averages of the six harmonic mode energies (calculated just as was done by 
Benettin^) are shown as functions of time for two different initial conditions, with total energies 
of 0.5 and 2.0. The six-particle chain of unit mass particles with least-energy coordinates of 
±0.5, ±1.5, and ± 2.5 is bounded by two additional fixed particles at ±3.5. 

Initially we choose the particles equally spaced and give all the energy to Particle 1, E — 
Pi/2. The left side of the Figure corresponds to an initial momentum of 1 while the right 
side follows a similarly long trajectory (100 million Runge-Kutta timesteps) starting with 
the initial momentum pi — 2. Fermi was surprised to find that at moderate energies there 
was no real tendency toward equilibration despite the anharmonic forces. Thus the averaging 
techniques of statistical mechanics can't usefully be applied to such oversimplified systems. 

Fermi also carried out some groundbreaking two-dimensional work. He solved Newton's 
equations of motion, 

{ mr = F(r) } , 

and didn't bother to describe the integration algorithm. A likely choice would be the time- 
reversible centered second- difference "Leapfrog" algorithm, 

{ r t+dt = 2r t - r t _ dt + {F/m) t {dtf } , 

where the timestep dt is a few percent of a typical vibrational period. The dominant error 
in this method is a "phase error", with the orbit completing prematurely. A harmonic 
oscillator with unit mass and force constant has a vibrational period of 2ir. The second- 
difference Leapfrog algorithm's period is 6, rather than 6.2832, for a relatively large timestep, 
dt = 1. A typical set of six (repeating) coordinate values for this timestep choice is: 

{ +2, +1,-1,-2,-1,+!, ... } 
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The Leapfrog algorithm diverges, with a period of 2\/2, as dt approaches y/2. 

Vineyard used the Leapfrog algorithm at the Brookhaven Laboratory, including irre- 
versible viscous quiet-boundary forces designed to minimize the effect of surface reflections 
on his simulations of radiation damage^. Alder and Wainwright, at the Livermore Labo- 
ratory, studied hard disks and spheres in parallel with Wood and Jacobsen's Monte Carlo 
work at the Los Alamos laboratory, finding a melting/freezing transition for spheres^ 1 ^. 
The disks and spheres required different techniques, with impulsive instantaneous momen- 
tum changes at discrete collision times. All these early simulations gave rise to a new 
discipline, "molecular dynamics" , which could be used to solve a wide variety of dynamical 
problems for gases, liquids, and solids, either at, or away from, equilibrium. By the late 
1960s the results of computer simulation supported a successful semiquantitative approach 
to the equilibrium thermodynamics of simple fluids- . 

In the 1970s Ashurst 17 (United States), Dremin^ (Union of Soviet Socialist Republics), 
Verletr^ (France), and Woodcock— (United Kingdom), were among those adapting molecular 
dynamics to the solution of nonequilibrium problems. Shockwaves, the subject of our third 
lecture, were among the first phenomena treated in the effort to understand the challenging 
problems of far-from-equilibrium many-body systems. 

C. Temperature Control a la Nose 

Shuichi Nose made a major advance in 1984 2 ->22, developing a dynamics, "Nose-Hoover 
dynamics", which provides sample isothermal configurations from Gibbs' and Boltzmann's 
canonical distribution, 

f(q,p) oc e~ n ^ kT ; H(q,p) = <%) + K{p) . 

The motion equations contain one or more friction coefficients {£} which influence the 
motion, forcing the longtime average of one or more of the pf to be mkTi. 

{ mfi =Fi- Qpi ; (i = [(p 2 i/mkTi) - l]/r 2 } . 

The thermostat variable ( can introduce or extract heat. The adjustable parameter r is 
the characteristic time governing the response of the thermostat variable (. A useful special 
case that follows from Nose's work in the limit r — > is "Gaussian" isokinetic dynamics, a 
dynamics with fixed, rather than fluctuating, kinetic energy K(p) = K , 
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In 1996 Dettmann showed that the Nose-Hoover equations of motion follow generally 
from a special Hamiltonian, without the need for the time scaling Nose used in his original 
work: 

H Dcttmann = s[$(q) + K{p/ a) + #kT\ns + #kT(p s r) 2 /2] = . 

Here the friction coefficient is ( = #kTr 2 p s , where p s is the Hamiltonian momentum con- 
jugate to s. The trick of setting the Hamiltonian equal to a special value, 0, is essential to 
Dettmann's derivation^. 

Consider the simplest interesting harmonic oscillator with unit mass, force con- 

stant, temperature, and relaxation time: 

U = s[q 2 + (p/s) 2 + In s 2 + p 2 s ]/2 = -> 

q= (p/s) ; p= -sq ; s = sp s ; p s = -[0] + {p/s) 2 - 1 -»- 

q = (l/s)p- (p/s)(s/s) = -q- (q ; p s = ( = q 2 - 1 . 

The time average of the ( equation shows that the longtime average of q 2 is unity. In 
particular applications r should be chosen to maximize the efficiency of the simulation by 
minimizing the necessary computer time. 

Runge-Kutta integration is a particularly convenient method for solving such sets of 
coupled first-order differential equations. The fourth-order method is the most useful. The 
time derivative is an average from four evaluations, { 2/0,2/1,2/2,2/3 }, of the righthand sides 
of all the differential equations, here collected in the form of a single differential equation 
for the vector y: 

2/i=2/o + (dt/2)y ; 2/2 = 2/0 + (dt/2)y 1 ; 2/3 = 2/o + dty 2 ; 

Vdt = 2/o + (dt/ 6) [y Q + 2y x + 2y 2 + 2/3] • 

The Runge-Kutta energy decays with time as dt 5 at a fixed time for a chosen timestep 
dt. Here the vector y is (q,p) so that 

y = (q,p) = (+p,-q) ■ 

For small dt the Runge-Kutta trajectory for a harmonic oscillator with the exact trajec- 
tory q = cos(t) has an error Sq = +dt t sin(£)/120. The corresponding Leapfrog error is 
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FIG. 2: Comparison of the maximum error (which occurs near a time of 3n/2), in a harmonic 
oscillator coordinate for the Leapfrog and Fourth-Order-Runge-Kutta integrators. The abscissa 
shows the logarithm of the number of force evaluations (which varies from about 20 to about 400) 
used during a full vibrational period, 2tt. The oscillator equations of motion are q = p ; q = ft = —q- 

5q = —dt 2 tsin(t)/2A. The two methods should give equally good solutions (where the two 
curves in the Figure cross) when 



corresponding to about 45 force evaluations per oscillator period . 

For a 14-digit-accurate trajectory calculation, with dtjut — 0.001 and dt^p = 0.00025, the 
Runge-Kutta error would be smaller than the Leapfrog error by seven orders of magnitude. 
At the cost of additional programming complexity choosing one of the fourth-order Gear 
integrators can reduce the integration error by an additional factor of ~ 60 25 . 

D. Connecting Microscopic Dynamics to Macroscopic Physics 

To connect the microscopic dynamics to macroscopic thermodynamics and continuum 
mechanics is quite easy for a homogeneous system confined to the volume V. A numerical 
solution of the equations of motion for the coordinates and momenta, { q,p }, makes it 
possible to compute the energy E, the temperature tensor T, the pressure tensor P, and the 
heat-flux vector Q: 




E = $(q) + K(p) = ^ <t>ij + E^/(2m) ; 



i<j i 



T xx = (p 2 x /mk) = ^2(p 2 Jmk)i/N ; T yy = (p 2 y /mk) ; 
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PV = F ij r ij + Y.(PP/ mk )i 5 

i<j i 

QV = J2 F ij ■ Vifij + J2(ep/mk)i . 

i<j i 

These expressions can be derived directly from the dynamics, by computing the mean mo- 
mentum and energy fluxes (flows per unit area and time) in the volume V. Alternatively 
they can be derived by multiplying the Newtonian equations of motion by (p/m) (giving 
the "Virial Theorem") or by e (giving the "Heat Theorem") and time averaging^. We will 
see that local versions of these definitions lead to practical implementations of numerical 
hydrodynamics at atomistic length and time scales. 

The thermomechanical bases of these relations are statistical mechanics and kinetic the- 
ory. Hamilton's mechanics yields Liouville's theorem for the time derivative of the many- 
body phase-space probability density following the motion: 

/ / / = din f I dt = [Hamiltonian Mechanics] ; 

Nose-Hoover mechanics opens up the possibility for / to change: 

f ' I f — dhaf /dt — C = —E/kT = S cxt /k [Nose — Hoover Mechanics] . 

The primary distinction between nonequilibrium and equilibrium systems lies in the friction 
coefficients { C }• At equilibrium (ordinary Newtonian or Hamiltonian dynamics) the average 
friction vanishes while in nonequilibrium steady states {J2 k() = S cxt > it is equal to the 
time-averaged entropy production rate. 

In any stationary nonequilibrium state the sum of the friction coefficients is necessarily 
positive - a negative sum would correspond to phase-space instability incompatible with 
a steady state. An important consequence of the positive friction is that the probability 
density for these states diverges as time goes on, indicating the collapse of the probability 
density onto a fractal strange attractor. Fractals differ from Gibbs' smooth distributions in 
that the density is singular, and varies as a fractional power of the coordinates and momenta 
m phase space^ 1 ^— . 

E. Fractal Phase-Space Distributions 

The harmonic oscillator problem is not ergodic with Nose-Hoover dynamics. One way 
to make it so is to fix the fourth moment of the velocity distribution as well as the second. 
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This improvement also makes it possible to study interesting nonequilibrium oscillator-based 
problems, such as the conduction of energy from hot to cold through the oscillator motion. 
Figure 3 shows the time development of (the two-dimensional projection of) such a problem. 
The isothermal oscillator, along with two friction coefficients, {C;£} ; fixing the second and 
fourth moments, {{p 2 ,p A )) has a Gaussian distribution in its four-dimensional phase space. A 
special nonisothermal case, with a coordinate-dependent temperature leading to heat flow, 
generates a 2.56-dimensional fractal in the four-dimensional {q,p, C 5 £} phase space. The 
dynamics governing this continuous nonequilibrium motion is as follows: 

q = p; p = -q - (p - £p 3 ; 

C = b 2 -T]; i=[p 4 -3p 2 T]; T = T(q) = 1 + tanh(g) . 

Here time averages of the control-variable equations show that the second and fourth mo- 
ments satisfy the usual thermometric definitions: 

(P 2 ) = (T) ; (p 4 ) = S(p 2 T) . 

The phase-space distribution for this oscillator has an interesting fractal nature^ 1 ^. Figure 
3 shows how the continuous trajectory comes to give a fractal distribution, as is typical 
of thermostated nonequilibrium problems. Besides the aesthetic interest that this model 
provides, it illustrates the possibilities for controlling moments of the velocity distribution 
beyond the first and second, as well as the possibility of introducing a coordinate-dependent 
temperature directly into the motion equations. 

Figure 4 shows a more typical textbook fractal, the Sierpinski sponge, in which the 
probability density is concentrated on a set of dimension 2.727. In almost all of the largest 
cube the density vanishes. Unlike the multifractal of Figure 3, the Sierpinski sponge is 
homogeneous, so that an ra-fold enlarged view of a small part of the sponge, with an overall 
volume 1/27™ of the total, looks precisely like the entire object. 

F. The Galton Board 

The situation with impulsive forces is quite different. Whenever impulsive collisions 
occur the phase-space trajectory makes a jump in momentum space, from one phase point 
to another. Consider the simplest interesting case: a single point mass, passing through a 
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FIG. 3: This (£, £) projection of the doubly-thermostated oscillator fractal is shown at five succes- 
sive stages of temporal resolution. The time intervals between successive points range from 0.001, 
the Runge-Kutta timestep, to 10.0, showing how a continuous trajectory can lead to a fractal 
object. 




FIG. 4: Sierpinski Sponge, constructed by removing 7 of the 27 equal cubes contained in the 
unit cube, leaving 20 smaller cubes, and then iterating this process ad infinitum leaving a 2.727- 
dimensional fractal of zero volume. 

triangular lattice of hard scatterers^^i. That model generates exactly the same ergodic 
dynamics as does a periodic two-hard-disk system with no center-of-mass motion: 

T\ + r 2 = ; Vi + v 2 = ■ 

By adding a constant field and an isokinetic thermostat to the field- dependent motion, the 
trajectory tends smoothly toward the field direction until a collisional jump occurs. Over 
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FIG. 5: A series of 200,000 Galton Board collisions are plotted as separate points, with ordinate 
— 1 < sin(/3) < 1 and abscissa < a < it, where a is measured relative to the field direction, as 
shown in Figure 6. 

long times (Figure 5 is based on 200,000 collisions) an extremely interesting nonequilibrium 
stationary state results, with a fractal phase-space distribution. The example shown in the 
Figure has an information dimension of 1.832. As a consequence, the coarse-grained entropy, 
— k(f In f), when evaluated with phase-space cells of size 5, diverges as S~ 0A68 , approaching 
minus infinity limiting case. 

The probability densities for nonequilibrium steady states, such as the Galton Board, 
shown in Figure 5, are qualitatively different to the sponge, where the probability density 
is equally singular wherever it is nonzero. The Galton Board's nonequilibrium probability 
density is nonzero for any configuration consistent with the initial conditions on the dy- 
namics. Further, the (multi)fractal dimension of these inhomogeneous distributions varies 
throughout the phase space. 

The concentrated nature of the nonequilibrium probability density shows first of all that 
nonequilibrium states are very rare in phase space. Finding one by accident has probability 
zero. The time reversibility of the equations of motion additionally shows that the probabil- 
ity density going forward in time contracts (onto a strange attractor), and so is necessarily 
stable relative to a hypothetical reversed trajectory going backward in time, which would 
expand in an unstable way. This symmetry breaking is a microscopic equivalent of the Sec- 
ond Law of Thermodynamics, a topic to which we'll return. It is evidently closely related to 
the many "fluctuation theorems"— 1 ^ which seek to give the relative probabilities of forward 
and backward nonequilibrium trajectories as calculated from Liouville's Theorem. 
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FIG. 6: The Galton Board geometry is shown, defining the angles a. and (3 identifying each 
collision. The unit cell shown here, extended periodically, is sufficient to describe the problem of a 
moving particle in an infinite lattice of scatterers. 

G. Determination of Transport Coefficients via NEMD 

With measurement comes the possibility of control. Feedback forces, based on the results 
of measurement, can be used to increase or decrease a "control variable" (such as the fric- 
tion coefficient ( which controls the kinetic temperature through a "thermostating" force). 
Equations of motion controlling the energy, or the temperature, or the pressure, or the heat 
flux, can all be developed in such a way that they are exactly consistent with Green and 
Kubo's perturbation-theory of transport^. That theory is a first-order perturbation theory 
of Gibbs' statistical mechanics. It expresses linear-response transport coefficients in terms 
of the decay of equilibrium correlation functions. For instance, the shear viscosity rj can be 
computed from the decay of the stress autocorrelation function: 



and the heat conductivity k can be computed from the decay of the heat flux autocorrelation 
function: 



Nose's ideas have made it possible to simulate and interpret a host of controlled nonequilib- 
rium situations. A Google search for "Nose-Hoover" in midJuly of 2010 produced over eight 
million separate hits. 

Figure 7 shows a relatively-simple way to obtain transport coefficients using nonequi- 
librium molecular dynamics. Ashurslr^, in his thesis work at the University of California, 
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FIG. 7: A Four Chamber viscous flow. Solid blocks (filled circles), move antisymmetrically to the 
left and right, so as to shear the two chambers containing Newtonian fluid (open circles). This 
geometry makes it possible to characterize the nonlinear differences among the diagonal components 
of the pressure and temperature tensors. 

"Dense Fluid Shear Viscosity and Thermal Conductivity via Molecular Dynamics", intro- 
duced two "fluid walls", with different specified velocities and/or temperatures, in order to 
simulate Newtonian viscosity and Fourier heat flow. Figure 7, a fully periodic variation of 
Ashurst's idea, shows two "reservoir" regions, actually "solid walls", separating two New- 
tonian regions. In both the Newtonian regions momentum and energy fluxes react to the 
different velocities and temperatures imposed in the "wall" reservoirs. This four-chamber 
technique produces two separate nonequilibrium profiles^— . 

In the Newtonian chambers, where no thermostat forces are exerted, the velocity or 
temperature gradients are nearly constant, so that accurate values of the viscosity and heat 
conductivity can be determined by measuring the (necessarily constant) shear stress or the 
heat flux: 

T) = -P xy /[{dv y /dx) + (dv x /dy)} ; k = -Q x /(dT/dx) . 
H. Nonlinear Transport 

This same "solid-wall" or "four-chamber" method has been used to study a more compli- 
cated aspect of nonequilibrium systems, the nonlinear contributions to the fluxes. Because 
the underlying phase-space distributions are necessarily fractal it is to be expected that 
there is no analytic expansion of the transport properties analogous to the virial (powers of 
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the density) expansion of the equilibrium pressure. Periodic shear flows, with the mean x 
velocity increasing linearly with y, 

{x = (p x /m) + ey ; y= (p y /m) } ; 

can be generated with any one member of the family of motion equations: 

{ fix = Fx - £Ot x p y - C,Px ; fiy = F y - ea y p x - (p y } , 

so long as the sum a x + a y is unity and ( is chosen to control the overall energy or temper- 
ature. Careful comparisons of the two limiting approaches, 

a x = ; a y = 1 [Doll's] ; 

a x — 1 ; a y = [s'lloD] , 

with corresponding boundary-driven four-chamber flows show that though both of the algo- 
rithms satisfy the nonequilibrium energy requirement: 

E = -eP xy V , 

exactly, neither of them provides the correct "normal stress" difference, P xx — P yy . 

This same problem highlights another interesting parallel feature of nonequilibrium sys- 
tems, the tensor nature of temperature^ - -^— . In a boundary-driven shearflow with the 
repulsive pair potential, 

0(r < 1) = 100(1 -r 2 ) 4 , 
the temperature tensors in the Newtonian regions show the orderings 

(pI) > (pI) > (pI) < — > Txx > T zz > T yy [Boundary Driven] . 

The homogeneous periodic shear flows generated with the Doll's and s'lloD algorithms show 
instead two other orderings: 

T xx > T yy > T zz [s'lloD] and T yy > T xx > T zz [Doll's] , 

so that neither the Doll's nor the s'lloD algorithm correctly accounts for the nonlinear prop- 
erties of stationary shear flows^. Nonequilibrium molecular dynamics provides an extremely 
versatile tool for determining nonlinear as will as linear transport. We will come back to 
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FIG. 8: Rayleigh-Benard problem, simulated with 5000 particles. The fluid-wall image particles 
which enforce the thermal and velocity boundary conditions are shown as circles above and below 
the main flow. 

tensor temperature in the third lecture, on Shockwaves. Nonlinear transport problems can 
require the definition of local hydrodynamic variables whenever the system is inhomoge- 
neous, as it is in boundary-driven shear and heat flows. 

Thermostats, ergostats, barostats, and many other kinds of constraints and controls 
simplify the treatment of complex failure problems with molecular dynamics. Using the 
Doll's and s'lloD ideas it is quite feasible to study the stationary nonequilibrium flow of 
solids, "plastic flow", in order to interpret nonsteady failure problems like fracture and 
indentation. Nonequilibrium molecular dynamics makes it possible to remove the irreversible 
heat generated by strongly nonequilibrium processes such as the machining of metals. The 
basic idea of control can be implemented from the standpoint of Gauss' Principle, which 
states that the smallest possible constraint force should be used to accomplish control. 
Near equilibrium a more reliable basis is Green and Kubo's linear-response theory. This 
can be used to formulate controls consistent with exact statistical mechanics in the linear 
regime, just as was done in deriving the Doll's and s'lloD approaches to simulating shear 
flow. 

A slightly more complex problem is illustrated in Figure 8. A nonequilibrium system 
with fixed mass is contained within two thermal "fluid wall" boundaries, hot on the bottom 
and cold on the top, with a gravitational field acting downward. If the gradients are small 
the fluid is stationary, and conducts heat according to Fourier's Law. When the Rayleigh 
Number, 

R = gL {dlnT /dy)/(vK) ; v = rj/p , 
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exceeds a critical value (which can be approximated by carrying out a linear stability analysis 
of the hydrodynamic equations) two rolls, one clockwise and the other counterclockwise 
provide another, faster, mode of heat transfer. At higher values of R the rolls oscillate 
vertically; at higher values still the rolls are replaced by chaotic heat plumes, which move 
horizontally. With several thousand particles molecular dynamics provides solutions in good 
agreement with the predictions of the Navier-Stokes-Fourier equations. 

This problem^— is specially interesting in that several topologically different solutions 
can exist for exactly the same applied boundary conditions. Carol will talk more about 
this problem in her exposition of Smooth Particle Applied Mechanics, "SPAM" . SPAM pro- 
vides a useful numerical technique for interpolating the particle properties of nonequilibrium 
molecular dynamics onto convenient spatial grids. 

II. PARTICLE-BASED CONTINUUM MECHANICS & SPAM 
A. Introduction and Goals 

Smooth Particle Applied Mechanics, "SPAM", was invented at Cambridge, somewhat 
independently, by Lucy and by Monaghan in 1977". The particles both men considered 
were astrophysical in size as their method was designed to treat clusters of stars. SPAM 
can be used on smaller scales too. SPAM provides a simple and versatile particle method 
for solving the continuum equations numerically with a twice-differentiable interpolation 
method for the various space-and-time-dependent field variables (density, velocity, energy, 
...) . SPAM looks very much like "Dissipative Particle Dynamics"-^, though, unlike DPD, it 
is typically fully deterministic, with no stochastic ingredients. Three pedagogical problems 
are discussed here using SPAM: the free expansion of a compressed fluid; the collapse of a 
water column under the influence of gravity; and thermally driven convection, the Rayleigh- 
Benard problem. Research areas well-suited to graduate research (tensile instability, angular 
momentum conservation, phase separation, and surface tension) are also described. 

SPAM provides an extremely simple particle-based solution method for solving the con- 
servation equations of continuum mechanics. For a system without external fields the basic 
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partial differential equations we aim to solve are: 

p = -pV • v ; 

pv = -V P ; 

pe = -Vv : P - V -Q . 

SPAM solves the equations by providing a particle interpretation for each of the continuum 
variables occuring in these conservation laws. The main difficulty in applying the method 
involves the choice and implementation of boundary conditions, which vary from problem 
to problem. 

B. SPAM Algorithms and the Continuity Equation 

The fluid dynamics notation here, {p, v, e, P, Q}, with each of these variables dependent 
on location r and time t, is standard but the SPAM particle interpretation of them is novel. 
The density p and momentum density (pv ) at any location r are local sums of nearby 
individual particle contributions, 

P( r ) = J2 m j w ( r - r j) ; P( r i) = ^mjw(ri - rj) ; p(r)v(r) = ^m^wfr - r,) , 

j j j 

where particles have an extent h, the "range" of the weight function w, so that only those 
particles within h of the location r contribute to the averages there. 

In the second expression (for the density at the particle location r^) the "self" term 
(ri = Tj) is included so that the two definitions coincide at the particle locations. The 
weight function w, which describes the spatial distribution of particle mass, or region of 
influence for particle j, is normalized, has a smooth maximum at the origin, and a finite 
range h, at which both w' and w" vanish. The simplest polynomial filling all these needs is 
Lucy's^, here normalized for two-dimensional calculations: 

w 2D {r <h) = (5/vr/i 2 )[l - 6a; 2 + 8a; 3 - 3x 4 ] ; x = r/h . 

Monaghan's weight function, shown for comparison in the Figure, uses two different poly- 
nomials in the region where w is nonzero. The range h of w(r < h) is typically a scalar, 
chosen so that a few dozen smooth particles contribute to the various field-point averages at 
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FIG. 9: Lucy's and Monaghan's weight functions. Both functions are normalized for two space 
dimensions and h = 3. The weight function w(r < h) describes the spatial influence of particles to 
properties in their neighborhood, as explained in the text. 

a point. As shown in Figure 9 Lucy's function looks much like a Gaussian, but vanishes very 
smoothly as r — > h. By systematically introducing the weight function into expressions for 
the instantaneous spatial averages of the density, velocity, energy, pressure, and heat flux, 
the continuum equations at the particle locations become ordinary differential equations 
much like those of molecular dynamics. The method has the desirable characteristic that 
the continuum variables have continuous first and second spatial derivatives. 

The continuity equation (conservation of mass) is satisfied automatically. At a fixed point 
r in space, the time derivative of the density depends upon the velocities of all those particles 
within the range h of r: 

(dp/dt) r = '^2'ITljVj ■ VjW r j = — ^ TTljVj • V ' r W r j , 
j 3 

where Vj is the velocity of particle j. On the other hand, the divergence of the quantity (pv ) 
at r is: 

j 

establishing the Eulerian and Lagrangian forms of the continuity equation: 

(dp/dt) r = -V r ■ (pv) < — > p = -pV • v . 

These fundamental identities linking the density and velocity definitions establish the 
smooth-particle method as the most "natural" for expressing continuous field variables in 
terms of particle properties. 
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FIG. 10: Contours of average density (middle row) and average temperature (bottom row) cal- 
culated from the instantaneous 16,384-particle snapshots (top row) taken during a free expansion 
simulation. The last picture in each row corresponds to two sound traversal times. 

The smooth-particle equations of motion have a form closely resembling the equations of 
motion for classical molecular dynamics: 

{ rrijVj = -J2 m j m k[( p /P 2 )j + ( P /P 2 )k\ ■ ^jW jk } • 

k 

It is noteworthy that the field velocity at the location of particle i 

v(r = r i )=Y / VjWij/ w ij = Yl m j v i w ij/p( r = r i) > 
j j j 

(where the "self" term is again included) is usually different to the particle velocity Vi, 
opening up the possibility for computing velocity fluctuations at a point, as we do in the 
next Section. 

Notice that the simple adiabatic equation of state P oc p 2 /2 gives exactly the same motion 
equations for SPAM as does molecular dynamics. That isomorphism pictures the weight 
function w(r) as the equivalent of a short-ranged purely-repulsive pair potential. Thus the 
continuum dynamics of a special two-dimensional fluid become identical to the molecular 
dynamics of a dense fluid with smooth short-ranged repulsive forces.- We consider this case 
further in applying SPAM to the free expansion problem in the next Section. 

C. Free Expansion Problem 

Figure 10 shows snapshots from a free expansion problem in which 16,384 particles, 
obeying the adiabatic equation of state P oc (p 2 /2), expand to fill a space four times that 
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FIG. 11: Equilibrated column for two system sizes. Five density contours are indicated by changes 
in plotting symbols. The arrows corresponding to the contours were calculated analytically from 
the continuum force-balance equation, dP/dy = —pg. 

of the initial compressed gas. This problem provides a resolution of Gibbs' Paradox (that 
the entropy increases by Nk\n4 while Gibbs' Liouville-based entropy, — k(lnf), remains 
unchanged)^ 1 ^. Detailed calculations show that the missing Liouville entropy is embodied 
in the kinetic-energy fluctuations. When these fluctuations are computed in a frame moving 
at the local average velocity, 

v(r) = J2 w rj v j/J2 w rj > 

j 'i 

the corresponding velocity fluctuations, ((v 2 ) — (v) 2 ) are just large enough to reproduce 
the thermal entropy. Most of the spatial equilibration occurs very quickly, in just a few 
sound traversal times. The contours of average density and average kinetic energy shown 
here illustrate another advantage of the SPAM averaging algorithm. The field variables are 
defined everywhere in the system, so that evaluating them on a regular grid, for plotting or 
analyses, is easy to do. 

These local velocity fluctuations begin to be important only when the adiabatic expansion 
stretches all the way across the periodic confining box so that rightward-moving fluid collides 
with its leftward-moving periodic image and vice versa. The thermodynamic irreversibility 
of that collision process, reproduced in the thermal entropy, is just sufficient for the reversible 
dynamics to reflect the irreversible entropy increase, NkhiA. A dense-fluid version of this 
dilute-gas free expansion problem appears in Bill's lecture on Shockwaves. 
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FIG. 12: Water Column collapse for three system sizes The computational time for this two- 
dimensional problem varies as the three-halves power of the number of particles used because 
corresponding times increase as yN while the number of interactions varies as N. 

D. Collapse of a Fluid Column 

Figure 11 shows the distribution of smooth particles in an equilibrated periodic water 
column in a gravitational fields. Figure 12 shows snapshots from the subsequent collapse of 
the water column when the vertical periodic boundaries are released. Both the equilibration 
shown in Figure 11 and the collapse shown in Figure 12 use the simple equation of state 
P = p 3 — p 2 , chosen to give zero pressure at unit density. Here the gravitational field 
strength has been chosen to give a maximum density of 2 at the reflecting lower boundary. 
Initially, the vertical boundaries are periodic, preventing horizontal motion. After a brief 
equilibration period, the SPAM density profile can be compared to its analytic analog, 
derived by integrating the static version of the equation of motion: 

dPj dy = -pg . 

The arrows in Figure 11, computed from the analytic static density profile, show excellent 
agreement with the numerical SPAM simulation. 

In smooth particle applied mechanics (SPAM) the boundary conditions are invariably 
the most difficult aspect of carrying out a simulations^. Here we have used a simple mirror 
boundary condition at the bottom of the column and a periodic boundary at the sides, in the 
vertical direction. When the vertical periodic boundary constraint is released, rarefaction 
waves create a tensile region inside the falling column. By varying the size of the smooth 
particles the resolution of the motion can be enhanced, as Figure 12 shows. With "mirror 
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FIG. 13: Instantaneous temperature (below) and density (above) contours for the two-roll 
Rayleigh-Benard problem. The stationary continuum solution (left) is compared to a SPAM snap- 
shot with 5000 particles (right). 

boundaries", elaborated in the next Section, more complicated situations can be treated. 
With mirrors there is an image particle across the boundary, opposite to each SPAM particle, 
with the mirror particle's velocity and temperature both chosen to satisfy the corresponding 
boundary conditions. 

E. Rayleigh-Benard Convection 

Figure 13 shows a typical snapshot for a slightly more complicated problem, the Rayleigh- 
Benard problem, the convective flow of a compressible fluid in a gravitational field with the 
temperature specified at both the bottom (hot) and top (cold) boundaries. The velocities at 
both these boundaries must vanish, and can be imposed by using mirror particles resembling 
the image charges of electricity and magnetism. A particularly interesting aspect of the 
Rayleigh-Benard problem is that multiple solutions of the continuum equations can coexist, 
for instance two rolls or four, with exactly the same boundary conditions^— Such work 
has been used to show that neither the entropy nor the entropy production rate allows one to 
choose "the solution" . Which solution is observed in practice can depend sensitively on the 
initial conditions. The Rayleigh-Benard problem illustrates the need for local hydrodynamic 
averages describing the anisotropics of two- and three-dimensional flows. 

SPAM provides an extremely useful interpolation method for generating twice- 
different iable averages from particle data. In the following lecture this method will be used 
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to analyze a dense-fluid molecular dynamics shockwave problem, where all of the thermo- 
mechanical variables make near-discontinuous changes linking an incoming cold state to an 
outgoing hot one. The continuous differentiable field variables provided by SPAM make it 
possible to analyze the relatively subtle nonlinear properties of such strongly nonequilibrium 
flow fields. 

SPAM is a particularly promising field for graduate research. In addition to the many 
possible treatments of boundaries (including boundaries between different phases), the con- 
servation of angular momentum (when shear stresses are present) and the tensile instability 
(where w acts as an attractive rather than repulsive force) and the treatment of surface 
tension all merit more investigation. For a summary of the current State of the Art see our 
recent book 6 . 

III. TENSOR-TEMPERATURE SHOCKWAVES VIA MOLECULAR DYNAMICS 
A. Introduction and Goals 

Shockwaves are an ideal nonlinear nonequilibrium application of molecular dynamics. 
The boundary conditions are purely equilibrium and the gradients are quite large. The 
shockwave process is a practical method for obtaining high-pressure thermodynamic data. 
There are some paradoxical aspects too. Just as in the free expansion problem, time- 
reversible motion, with constant Gibbs' entropy, describes a macroscopically irreversible 
process in which entropy increases. The increase is third-order in the compression, for 
weak shocks 50 . The shockwave problem is a compelling example of Loschmidt's reversibility 
paradox. 

We touch on all these aspects of the shockwave problem here. We generate and analyze 
the pair of Shockwaves which results from the collision of two stress-free blocks^. The 
blocks are given initial velocities just sufficient to compress the two cold blocks to a hot 
one, at twice the initial density. Further evolution of this atomistic system, with the initial 
kinetic energy of the blocks converted to internal energy, leads to a dense-fluid version of 
the free expansion problem discussed earlier for an adiabatic gas. Here we emphasize the 
dynamical reversibility and mechanical instability of this system, show the shortcomings of 
the usual Navier-Stokes- Fourier description of Shockwaves, and introduce a two-temperature 
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continuum model which describes the strong shockwave process quite well. 
B. Shockwave Geometry 

There is an excellent treatment of Shockwaves in Chapter IX of Landau and Lifshitz' 
"Fluid Mechanics" text^. A stationary shockwave, with steady flow in the x direction, 
obeys three equations for the fluxes of mass, momentum, and energy derived from the three 
continuum equations expressing the conservation of mass, momentum, and energy: 

pv = p c u s = p H {u s - u p ) ; 

P xx + pv 2 = P c + Pcu 2 s = Ph + Ph(u s - u p ) 2 ; 

pv[e + (P xx /p) + (v 2 /2)} + Q x = 

[e + (P xx /p)}c + (u 2 s /2) = [e + {P xx /p)] H + (it, - u p ) 2 /2 . 

Figure 14 illustrates the shockwave geometry in a special coordinate frame. In this frame 
the shockwave is stationary. Cold material enters from the left at the "shock speed" u s 
and hot material exits at the right, at speed u s — u p , where u p is the "particle" or "piston" 
velocity. The terminology comes from an alternative coordinate system, in which motionless 
cold material is compressed by a piston (moving at u p ), launching a shockwave (moving at 

U a ). 

Eliminating the two speeds from the three conservation equations gives the Hugoniot 
equation, 

e H - e c = {Ph + Pc)(V c ~ V H )/2 , 

which relates the equilibrium pressures, volumes, and energies of the cold and hot states. 
Evidently purely equilibrium thermodynamic equation of state information can be obtained 
by applying the conservations laws to optical or electrical velocity measurements in the 
highly-nonequilibrium shockwave compression process. Ragan described the threefold com- 
pression of a variety of materials (using an atomic bomb explosion to provide the pressure) 
at pressures up to 60 Megabars, about 15 times the pressure at the center of the earth^. 

Hoover carried out simulations of the shockwave compression process for a repulsive 
potential, <f>(r) = r -12 , in 196 7 52 , but put off completing the project for several years, 
until computer storage capacity and execution speeds allowed for more accurate work—. 
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FIG. 14: Stationary shockwave in the comoving frame. Cold material enters at the left, with 
velocity +u s , and is decelerated by the denser hotter material which exits at the right, with 
velocity u s — u p . It is in this coordinate frame that the fluxes given in the text are constant. 
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FIG. 15: A series of snapshots showing the stability of a planar shockwave. Note that the decay 
of the initial sinewave profile is slightly underdamped. Here Dt = 2000<it is the time required for a 
shockwave to traverse the width shown here, 2000 Runge-Kutta timesteps with dt = 0.02/u s ~ 0.01. 

Comparison of Klimenko and Dremin's computer simulations^ using the Lennard- Jones 
potential, <f>(r) = r~ 12 — 2r~ 6 , showed that relatively weak Shockwaves (30 kilobars for 
argon, 1.5-fold compression) could be described quite well^ 1 ^ with the three-dimensional 
Navier-Stokes equations, using Newtonian viscosity and Fourier heat conduction: 

P = P eg - AV ■ v - r][Vv + W] ; 

A 2 £> = Vv - V ; A 3D = t] V - (2/3)77 ! 
Q = - K VT . 

Here A is the "second viscosity" , defined in such a way that the excess hydrostatic pressure 
due to a finite strain rate is — r]y V • v. The shear viscosity rj and heat conductivity k were 
determined independently using molecular dynamics simulations. The small scale of the 
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FIG. 16: Snapshots near the beginning (upper) and end (lower) of an inelastic collision between 
two 400-particle blocks. The initial velocities, ±0.965, are just sufficient for a twofold compression 
of the cold material. The unit-mass particles interact with a short-ranged potential (10/-7r)(l — r) 3 
and have an initial density 



waves^, just a few atomic diameters, was welcomed by high-pressure experimentalists weary 
of arguing that their explosively-generated Shockwaves measured equilibrium properties. 

It is necessary to verify the one-dimensional nature of the waves too. It turns out that 
Shockwaves do become planar very rapidly, at nearly the sound velocity. The rate at which 
sinusoidal perturbations are damped out has been used to determine the plastic viscosity of 
a variety of metals at high pressure^. Figure 15 shows the rapid approach to planarity of a 
dense-fluid shockwave 9 . 

Stronger Shockwaves, where the bulk viscosity is more important (400 kilobars for ar- 
gon, twofold compression), showed that the Navier-Stokes description needs improvement 
at higher pressures. In particular, within strong Shockwaves temperature becomes a sym- 
metric tensor, with T xx » T yy , where x is again the propagation direction. In addition, 
the Navier-Stokes-Fourier shockwidth, using linear transport coefficients, is too narrow. The 
tensor character of temperature in dilute-gas Shockwaves had been carefully discussed in the 
1950s by Mott-Smith^.. 

C. Analysis of Instantaneous Shockwave Profiles using SPAM Averaging 

Data for systems with impulsive forces, like hard spheres, require both time and 
space averaging for a comparison with traditional continuum mechanics. Analy- 
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ses of molecular dynamics data with continuous potentials need no time averag- 
ing, but still require a spatial smoothing operation to convert instantaneous particle 
data, {x,y,p x ,p y }i, including {P, Q,T, e}j, to equivalent continuous continuum profiles, 
{p(r, t), v{r, t), e(r, t), P{r, t),T(r, t), Q{r, t)}. 

The potential parts of the virial-theorem and heat-theorem expressions for the pressure 
tensor P and the heat-flux vector Q, 

PV = /' + Y(pp/mk)i ; 

i<j i 

QV = ^ F ij ■ Vifij + J2(ep/mk)i , 

i<j i 

can be apportioned in at least three "natural" ways between pairs of interacting particles^* 6 -. 

Consider the potential energy of two particles, |ti2 | ) - This contribution to the system's 
energy can be split equally between the two particle locations, r± and r 2 , or located at the 
midpoint between them, (r\ + r 2 )/2, or distributed uniformly 56 along the line r\ — r 2 joining 
them. These three possibilities can be augmented considerably in systems with manybody 
forces between particles of different masses. It is fortunate that for the short-ranged forces we 
study here the differences among the three simpler approaches are numerically insignificant. 
Once a choice has been made, so as to define particle pressures and heat fluxes, these can 
in turn be used to define the corresponding continuum field variables at any location r by 
using the weight-function approach of smooth particle applied mechanics: 

p ( r ) =J2 p j w rj/Y,w rj ; Q(r) = J2QjW rj /J2w rj . 

j j j j 

By using this approach our own simulations have characterized another constitutive com- 
plication of dense-fluid Shockwaves - the time delays between [1] the maximum shear stress 
and the maximum strainrate and [2] the maximum heat flux and and the maxima of the 
two temperature gradients (dT xx /dx) and (dT yy / 'dx)^-. The study of such delays goes 
back to Maxwell. The "Maxwell relaxation" of a viscoelastic fluid can be described by the 
modelL&JiL 

a + TO = T]€ . 

so that stress reacts to a changing strainrate after a time of order r. Cattaneo considered the 
same effect for the propagation of heat. The phenomenological delays, found in the dynam- 
ical results, are a reminder that the irreversible nature of fluid mechanics is fundamentally 
different to the purely- reversible dynamics underlying it. 
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FIG. 17: Shock Thermal and Mechanical Profiles from molecular dynamics are shown at the top. 
Corresponding numerical solutions of the generalized continuum equations are shown at the bottom. 
This rough comparison suggests that the generalized equations can be fitted to particle simulations. 
The generalized equations use tensor temperature and apportion heat and work between the two 
temperatures T xx and T yy . They also include delay times for shear stress, for heat flux, and for 
thermal equilibration. 

The irreversible shock process is particularly interesting from the pedagogical standpoint. 
The increase in entropy stems from the conversion of the fluid's kinetic energy density, pv 2 /2 
to heat. To avoid the need for discussing the work done by moving pistons of Figure 14, 
we choose here to investigate Shockwaves generated by symmetric collisions of two stressfree 
blocks, periodic in the direction parallel to the shockfront. The entropy increase is large 
here (a zero-temperature classical system has an entropy of minus infinity). Figure 16 shows 
two snapshots for a strong shockwave yielding twofold compression of the initial cold zero- 
pressure lattice. The mechanical and thermal variables in a strong dense-fluid shockwave 
are shown in Figure 17. In order to model these results two generalizations of traditional 
hydrodynamics need to be made: the tensor nature of temperature and the delayed response 
of stress and heat flux both need to be treated. A successful approach is described next. 

D. Macroscopic Generalizations of the Navier-Stokes-Fourier Approach 

By generalizing continuum mechanics to include tensor temperature and the time delays 
for stress and heat flux, 

a + T& = i]i ; Q + tQ = -kVT , 
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FIG. 18: Shear Stress lags behind the strainrate. The molecular dynamics gradients, using smooth- 
particle interpolation, are much more sensitive to the range of the weighting function than are the 
fluxes. The results here are shown for h = 2, 3, 4, with line widths corresponding to h. 
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FIG. 19: Heat Flux lags behind the temperature gradients. The molecular dynamics gradients, 
using smooth-particle interpolation, are much more sensitive to the range of the weighting function 
than are the fluxes. The results here are shown for h = 2, 3, 4, with line widths corresponding to 
h. 

with an additional relaxation time describing the joint thermal equilibration of T xx and 
T yy to a common temperature Th the continuum and dynamical results can be made 
consistent£~-i2£>2IL^. In doing this we partition the work done and the heat gained into 
separate longitudinal (x) and transverse (y) parts: 

pf xx oc -aVt; : P - /3V ■ Q + piTyy - T xx )/r ; 

pfyy cx -(1 - a)Vv : P - (1 - /3)V • Q + p{T xx - T yy )/r . 
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Solving the time-dependent continuum equations for such shockwave problems is not 
difficult^ 1 ^. If all the spatial derivatives in the continuum equations are expressed as centered 
differences, with density defined in the center of a grid of cells, and all the other variables 
(velocity, energy, stress, heat flux, ...) at the nodes defining the cell vertices, fourth-order 
Runge-Kutta integration converges nicely to solutions of the kind shown in Figure 17. 

E. Shockwaves from Two Colliding Blocks are Nearly Reversible 

To highlight the reversibility of the irreversible shockwave process let us consider the col- 
lision of two blocks of two-dimensional zero-pressure material, at a density of y^4/3 (nearest- 
neighbor distance is unity, as is also the particle mass). Measurement of the equation of 
state with ordinary Newtonian mechanics, using the pair potential, 

<f>{r < 1) = (10/tt)(1 -rf , 

indicates (and simulation confirms) that the two velocities u s and u p , 

u s = 2u p = 1.930 , 

correspond to twofold compression with a density change y4/3 — > 2^4/3. To introduce a 
little chaos into the initial conditions random initial velocities, corresponding to a temper- 
ature 10~ 10 were chosen. Because the initial pressure is zero the conservation relations are 
as follows: 

pv = v /4/3 x 1.930 = 2.229 ; 

Pxx + pv 2 = \A73 x 1.930 2 = 4.301 ; 

pv[e + (P xx /p) + pv 2 /2)] + Q X = ^4/3 x 1.930 3 /2 = 4.151 . 

Although the reversibility of the dynamics cannot be perfect, the shockwave propagates 
so rapidly that a visual inspection of the reversed dynamics shows no discrepancies over 
thousands of Runge-Kutta timesteps. To assess the mechanical instability of the shock 
compression process we explore the effects of small perturbations to the reversible dynamics 
in the following Sections. We begin by illustrating phase-space instability^^ for a simpler 
problem, the harmonic chain. 
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F. Linear Growth Rates for a Harmonic Chain 



Even the one-dimensional harmonic chain, though not chaotic, exhibits linear phase- 
volume growth in certain phase-space directions. Consider the equations of motion for a 
periodic chain incorporating an arbitrary scalefactor s +2 : 

{q = ps + 2 ; p = ( q+ - 2q + qJ)s~ 2 } ; 

the subscripts indicate nearest-neighbor particles to the left and right. The motion equations 
for a 2iV-dimensional perturbation vector S = (Sq, Sp) follow by differentiation: 

{ Sq = Sps +2 ; dp = (5q + - 25q + Sq-)s~ 2 } . 

If we choose the length of the perturbation vector equal to unity, the logarithmic growth 
rate, A = (d\nS/dt) qiP , is a sum of the individual particle contributions: 

A(<y) = Y.1 6( i 6 p( s+2 - 2s ~ 2 ) + Wq+ + Sq-)s- 2 } . 

For a large scale factor s +2 it is evident that choosing equal components of the vector provides 
the maximum growth rate, 



{Sq = 5p= Vl/2iV} -> A max = 2-V 2 . 
For s 2 small, rather than large, alternating signs give the largest growth rate, with 

{ +5q even = +5p odd = -5q odd = Sp cvcn } , 

the growth rate is 

The growth rate is 2~ 1//2 at the transition between the two regions, where s 2 > 2 1 / 2 . 

These same growth-rate results can be found numerically by applying "singular value 
decomposition" to the dynamical matrix D^^-. This analysis details the deformation of an 
infinitesimal phase-space hypersphere for a short time dt. During this time the hypersphere 
has its components Sq, 5p changed by the equations of motion: 

5 A (/ + Ddt) ■ 5 , 
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so that the growth and decay rates can be found from the diagonal elements of the singular 
value decomposition 



I + Ddt = U -W -V* ->{A = (l/dt)\nW) . 



Numerical evaluation gives the complete spectrum of the growth and decay rates. The 
maximum matches the analytic results given above. Although locally the growth rates 
{A(r, t)} are nonzero, the harmonic chain is not at all chaotic and the long-time-averaged 
Lyapunov exponents {A = (A(r, t))}, all vanish. Let us now apply the concepts of phase- 
space growth rates {A} and the Lyapunov exponents {A} to the shockwave problem. 

G. Linear Instability in Many Body Systems, A for Shockwaves 

The time reversibility of the Hamiltonian equations of motion guarantees that any station- 
ary situation shows both a long-time-averaged and a local symmetry between the forward 
and reversed directions of time. In such a case the N nonzero time-averaged Lyapunov 
exponents as well as the local growth rates, obey the relations 



The instantaneous Lyapunov exponents {A(t)} depend on the dynamical history, while the 
instantaneous diagonalized phase-space growth rates, which we indicate with A(t) rather 
than X(t), do not. 

The rates {A(t)} for different directions in phase space can be calculated efficiently from 
the dynamical matrix D, by using singular value decomposition, just as we did for the 
harmonic chain: 



Here we analyze a 480-particle shockwave problem, the collision of two blocks with x 
velocity components ±0.965. Figure 20 shows those particles making above-average con- 
tributions to the maximum phase-space growth rate at times 2, 4, and 6. At time 12 the 
particle velocities are all reversed, so that the configurations at times 22, 20 and 18 closely 
match those at times 2, 4, and 6. Generally there is six-figure agreement between the coor- 
dinates going forward in time and those in the reversed trajectory at corresponding times. 



{ X^+i-k + Afc} — ; { A N+ i_ k + A k } — . 



dq/dq 
^ dp/dq 
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FIG. 20: Phase-space Growth Rates {A} during the collision of two 240-particle blocks of length 
20. The collision leads to twofold compression of the original cold material at a time of order 
20/1.93 ~ 10.4. At time = 12 the velocities were reversed, so that the configurations at times 2, 
4, and 6 correspond closely to those at 22, 20, and 18 respectively. Those particles making above 
average contributions to the largest phase-space growth direction are indicated with open circles. 
The fourth-order Runge-Kutta timestep is dt = 0.002. 

Note this symmetry in Figure 20, where the most sensitive particles going forward and back- 
ward are exactly the same at corresponding times. The forward-backward agreement could 
be made perfect by following Levesque and Verlet's suggestion^ to use integer arithmetic in 
evaluating a time-reversible (even bit-reversible^ algorithm such as 

lnt[q t+dt - 2q t + q t ^ dt ] = lnt[F t dt 2 / m] or 

lnt[q t+ 2dt ~ qt+dt ~ qt-dt + qt-2dt] = lnt[(dt 2 / 'Am)(hF t+dt + 2F t + 5F t -dt)] . 

Evidently, as would be expected, from their definition, the point-function growth rates 
{A(r(i))} can show no "arrow of time" distinguishing the backward trajectory from the 
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forward one. We turn next to the Lyapunov exponents, which can and do show such an 
arrow. 

H. Lyapunov Spectrum in a Strong Shockwave 

Most manybody dynamics is Lyapunov unstable, in the sense that the length of the phase- 
space vector joining two nearby trajectories has a tendency to grow at a (time-dependent) 
rate \i(t) (with the time-averaged result Ai = (Ai(t)) > 0). Likewise, the area of a moving 
phase-space triangle, with its vertices at three nearby trajectories, grows at \i(t) + A 2 (£), 
with a time-averaged rate Ai + A 2 . The volume of a tetrahedron defined by four trajectories 
grows as Ai(t) + A 2 (t) + Xz(t), and so on. By changing the scale factor linking coordinates 
to momenta - the s +2 of the last Section - these exponents can be determined separately in 
either coordinate or momentum space. 

Posch and Hoover, and independently Goldhirsch, Sulem, and Orszag, discovered a 
thought-provoking representation of local Lyapunov exponents^^i. If an array of Lagrange 
multipliers is chosen to propagate a comoving corotating orthonormal set of basis vectors 
centered on a phase space trajectory, the diagonal elements express local growth and de- 
cay rates. These are typically quite different (and unrelated) in the forward and backward 
directions of time. 

Let us apply the Lyapunov spectrum^— to the phase-space instability of a strong shock- 
wave. Because the Lyapunov exponents, {A(t)}, are evaluated so as to reflect only the past, 
times less than t, we expect to find that the Lyapunov vector corresponding to maximum 
growth soon becomes localized near the shock front. Starting out with randomly oriented 
vectors the time required for this localization is about 1/2. The time-linked disparities be- 
tween the forward and backward motions suggest that the Lyapunov exponents can provide 
an "Arrow of Time" because the stability properties forward in time differ from those in the 
backward (reversed) direction of time^. 

Figure 21 shows the particles making above-average contributions to the largest of the 
local Lyapunov exponents, Ai(i). There are many more of these particles than the few which 
contribute to the largest of the phase-space growth rates, Ai(t). The shockwave simulation 
was run forward in time for 6000 timesteps, after which the velocities were reversed. The 
phase-space offset vectors, chosen randomly at time and again at time 12, became localized 
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time = 2 
time = 4 
time = 6 
time = 18 
time = 20 
time = 22 

FIG. 21: Phase-space Growth Rates {A} during the collision of two similar blocks which lead to 
the twofold compression of the original cold material. At time = 12 the velocities were reversed, 
so that the configurations at times 2 and 4 correspond to those at 22 and 20, respectively. The 
particles making above average contributions to the largest Lyapunov exponent are indicated with 
open circles. 

near the shockfront at a time of order 0.5. The particles to which the motion is most sensitive, 
as described by the Lyapunov exponent Xi(t) are more localized in space in the forward 
direction of time than in the backward direction. Evidently the Lyapunov vectors are more 
useful than the vectors corresponding to local growth rates in describing the irreversibility 
of Hamiltonian systems. 

IV. CONCLUSION 

Particle dynamics, both NEMD and SPAM, provides a flexible approach to the simulation, 
representation, and analysis of nonequilibrium problems. The two particle methods are 
closely related, making it possible to infer constitutive relations directly from atomistic 
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simulations. These useful tools provide opportunities for steady progress in understanding 
far-from-equilibrium states. It is our hope that these tools will become widely adopted. 
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